package scu.maqiang.tiangong;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3T41;
import scu.maqiang.fes.FES3T42;
import scu.maqiang.mesh.Mesh3T4;
import scu.maqiang.mesh.MeshFormat;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class MainEngineHeating {

	public static void main(String[] args) {
		String dir = "TianGong\\MainEngine\\ComputeData\\";
		String meshFile = dir + "MainEngine.mesh";
		Mesh3T4 mesh = new Mesh3T4().readMesh(meshFile, MeshFormat.mesh);
		ScalarFunc fixedRegion = (xyz, label, param) -> xyz[0] < 0.000001?1.0:-1.0;
		int fixedLabel = 10;
		mesh.setBoundaryLabel(fixedRegion, null, fixedLabel);
		
		FES3T41 fsHeat = new FES3T41(mesh);
		int dofHeat = fsHeat.getNdof();
		double kxx = 217.0;
		double rho = 2700;
		double specificHeat = 800;
		double theta = 0.5;
		double endtime = 100;
		double dt = 2.0;
		SRMatrix A = new SRMatrix(dofHeat);
		SRMatrix M = new SRMatrix(dofHeat);
	    fsHeat.assembleStiff(new double[] {kxx}, BVPType.COMMON, A);
	    fsHeat.assembleMass(new double[] {rho * specificHeat}, BVPType.CONSIST_MASS, M);
	    
	    SRMatrix LMatrix = M.add(theta * dt, A);
		SRMatrix RMatrix = M.add((theta - 1) * dt, A); 
		
		double[] heatFlux0 = MVO.readArray(dir + "MainEngineHeatflux0.dat");
		double[] heatFlux10 = MVO.readArray(dir + "MainEngineHeatflux10.dat");
		double[] heatFluxm10 = MVO.readArray(dir + "MainEngineHeatflux-10.dat");
		double[] heatFlux20 = MVO.readArray(dir + "MainEngineHeatflux20.dat");
		double[] heatFluxm20 = MVO.readArray(dir + "MainEngineHeatflux-20.dat");
		
		double[] pressure0 = MVO.readArray(dir + "MainEnginePressure0.dat");
		double[] pressure10 = MVO.readArray(dir + "MainEnginePressure10.dat");
		double[] pressurem10 = MVO.readArray(dir + "MainEnginePressure-10.dat");
		double[] pressure20 = MVO.readArray(dir + "MainEnginePressure20.dat");
		double[] pressurem20 = MVO.readArray(dir + "MainEnginePressure-20.dat");
		
		double E = 71.8e9;
		double Nu = 0.33;
		double alpha = 23.0e-6;
		double beta = E / (1 - 2 * Nu) * alpha;
		FES3T42 fsElastic = new FES3T42(mesh);
		int dofElastic = fsElastic.GetNdof();
		SRMatrix K = new SRMatrix(dofElastic);
		fsElastic.assembleStiff(new double[] {E, Nu}, BVPType.COMMON, K);
		
		double[] x = new double[dofHeat];
		double[] FFlux0 = new double[dofHeat];
		double time = 0.0;
		String resultFile = dir + "MainEngineHeating10.dat";
		mesh.toTecplot(resultFile, time, x);
		double[][] vecCoef = new double[1][];
		vecCoef[0] = heatFlux10;
		fsHeat.assembleFlux(vecCoef, BVPType.COMMON, FFlux0, 0, 10);
		Arrays.setAll(FFlux0, i -> FFlux0[i] * dt);
		double[] Rx = new double[dofHeat];
		double[] RHS = new double[dofHeat];
		NewIterSSolver heatSolver = new NewIterSSolver(LMatrix);
		System.out.println("Start Computing!");
		
		fsElastic.applyBC_MBN(K, Direct.All, 10);
		NewIterSSolver elasticSolver = new NewIterSSolver(K);
		//DirectSSolver elasticSolver = new DirectSSolver(K, fsElastic.getnBand());
		double[] RHS_elastic = new double[dofElastic];
		double[] RHS_Pressure = new double[dofElastic];
		vecCoef[0] = pressure10;
		fsElastic.assemblePressure(vecCoef, BVPType.COMMON, RHS_Pressure, 0, 10);
		vecCoef[0] = x;
		double[] d = new double[dofElastic];
		String resultFile2 = dir + "MainEngineDeforming10.dat";
		double[][] disp = new double[3][dofHeat];
		mesh.toTecplot(resultFile2, time, disp);
		String resultFile21 = dir + "MainEngineDeforming10_1.dat";
		mesh.toTecplotDeform(resultFile21, time, disp, 10.0);
		double[][] strain = new double[6][dofHeat];
		double[][] stress = new double[6][dofHeat];
		double[] misesStress = new double[dofHeat];
		String resultFile3 = dir + "MainEngineStain10.dat";
		String resultFile4 = dir + "MainEngineStress10.dat";
		double[][] stressOutput = new double[][]{ stress[0], stress[1], stress[2], stress[3], stress[4], stress[5], misesStress};
		mesh.toTecplot(resultFile3, time, strain);
		mesh.toTecplot(resultFile4, time, stressOutput);
		while(time < endtime) {
			time += dt;
			System.out.println("time = " + time);
			RMatrix.mul(x, Rx);
			Arrays.setAll(RHS, i -> Rx[i] + dt * FFlux0[i]);
			heatSolver.PCGSSOR(RHS, x, 1.65, 1);
			mesh.toTecplot(resultFile, time, x);
			
			fsElastic.assembleHeatLoad(new double[] {beta}, vecCoef, BVPType.COMMON, RHS_elastic);
			Arrays.setAll(RHS_elastic, i -> RHS_elastic[i] + RHS_Pressure[i]);
			//System.out.println(MVO.max(RHS_elastic) +"\t " + MVO.min(RHS_elastic));
			//System.exit(0);
			fsElastic.applyBC_MBN(RHS_elastic, Direct.All, 0.0, 10);
			elasticSolver.PCGSSOR(RHS_elastic, d, 1.5, 1);
			fsElastic.extractComponent(d, disp);
			mesh.toTecplot(resultFile2, time, disp);
			mesh.toTecplotDeform(resultFile21, time, disp, 10.0);
			Arrays.fill(RHS_elastic, 0.0);
			
			fsElastic.computeStrainStressEqStress(BVPType.COMMON, new double[] {E, Nu},  d, strain, stress, misesStress);
			mesh.toTecplot(resultFile3, time, strain);
			mesh.toTecplot(resultFile4, time, stressOutput);
			System.out.println("Strain: ");
			System.out.println("ex min: " + MVO.min(strain[0])+ " \tmax: " + MVO.max(strain[0]));
			System.out.println("ey min: " + MVO.min(strain[1])+ " \tmax: " + MVO.max(strain[1]));
			System.out.println("ez min: " + MVO.min(strain[2])+ " \tmax: " + MVO.max(strain[2]));
			System.out.println("exy min: " + MVO.min(strain[3])+ " \tmax: " + MVO.max(strain[3]));
			System.out.println("eyx min: " + MVO.min(strain[4])+ " \tmax: " + MVO.max(strain[4]));
			System.out.println("exz min: " + MVO.min(strain[5])+ " \tmax: " + MVO.max(strain[5]));
			
			System.out.println("Stress: ");
			System.out.println("tx min: " + MVO.min(stress[0])+ " \tmax: " + MVO.max(stress[0]));
			System.out.println("ty min: " + MVO.min(stress[1])+ " \tmax: " + MVO.max(stress[1]));
			System.out.println("tz min: " + MVO.min(stress[2])+ " \tmax: " + MVO.max(stress[2]));
			System.out.println("txy min: " + MVO.min(stress[3])+ " \tmax: " + MVO.max(stress[3]));
			System.out.println("tyx min: " + MVO.min(stress[4])+ " \tmax: " + MVO.max(stress[4]));
			System.out.println("txz min: " + MVO.min(stress[5])+ " \tmax: " + MVO.max(stress[5]));
			
			System.out.println("Von Mises Stress: ");
			System.out.println("min: " + MVO.min(misesStress)+ " \tmax: " + MVO.max(misesStress));
		}
	}

}
